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Preface 


The work described in this report was performed by the Information 
Systems Division of the Jet Propulsion Laboratory. 
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Abstract 

Modified divided differences (MDD) provide a good way of representing 
a polynomial passing through points with unequally spaced abcissas. 
This report gives recurrence relations for computing coefficients in 
either the monomial or Chebyshev basis from the MDD coefficients, and 
for computing the MDD coefficients for either the differentiated or 
the integrated polynomial. The latter operation is likely to be 
useful if MDD are used in a method for solving stiff differential 
equations. 
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1. Introduction 


A modification of divided differences was first suggested by 
Blanch [l] as a means of getting some of the desirable characteristics 
of differences when working with unequal intervals. Modified divided 
differences of the form used here were introduced in [2] and developed 
further in [s] . Shampine and Gordon use them in [4] and Jackson (5] 
has made a careful study of related computational approaches, concluding 
that the form used here is best. Jackson has also showed that the 
round-off characteristics of these modified divided differences are 
excellent. 


The polynomial of lowest degree passing through the points 
(t^, w(t^)), i = n, n-1, ..., n-q+1 can be expressed in terms of the 
MDD coefficients (n) according to the formula 


( 1 ) 




q-1 

S 4>.{n) c. (t) 
i=0 ’ 


Where we use notation similar to that in [3]. 
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( 2 ) 


^ - ^-1 


T - (t-t,)/h„ 


(>i) 1 + . . . + h„ . = t„-t„ . 1 = C,' i(n-l) + 

^ 1 ' ' n n-1 n-i n n-i-1 ^i-r ' n 

ai(n) = h^/?.(n) 

eo(f’) = 1 

6,-+l(n) = {Ci*(n)/Ci(n-l))B^.(n) 


n) = w(t^) 


(|>i+i(n) = 4>j(n) - B^{n)4).{n-1) 

‘^0,n " ^ 

‘^l.n = = 

■^l.r, “ t C,.2(n)/C^.,(n)] 1 2 2 

Here we think of the 4>^(n) as coefficients of a polynomial expressed in terms 
of the basis polynomials c. „(t). In [3], the orientation was focused more 
on the 4>^{n) which are the modified divided differences of w computed at t=t^. 


2. Converting to Coefficients for Other Basis Polynomials 

Salzer [6] gives a general recurrence for converting coefficients for 
one set of polynwial basis functions to coefficients for another set of poly- 
nomial basis functions. In Salzer's notation the c^.(t) (Henceforth we no 
longer explicitly include n when referring to the various parameters defined 
in eq. (?)) satisfy the recurrence. 


( 3 ) 


c^i(t) + (a(m) + b(m)x) c^(x) + c(m) c^_i(t) = 0 


where a(0) = 0, a(m) = ni > 1 


b(m) = -a, 
c(m) = 0 


m 
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i' 


1 


and for the monomial basis, Q,^(t) = t"’, A(m) = C(m) = 0, B{m) = -1. Simply 
substituting in his eq. (22) yields (his n = q-1) 

° *q-l 

4^^ ° Vk-2 * 4''* 

(k+1) (k) (k) 

\ ' “q-k-2'm-l '^^S-k-3^S-k-2^ % , m = 1. 2, k 

(k+1) (k) 

®k+l °‘q-k-2 ^k 


where k = 0, 1, ...» q-2, and C_-| = 0. The desired coefficients are then 
a^^ \ m=0,l,..., q-1. Since Uq = 1 , m > 0. Substitu- 
ting (^,pi|/aq.|^. 2 ) simplifies to 


(5) a^ - a^.24>q.-| 


(k+1) 


(k) 

*1 

°‘q-k-3*^q-k-2 

^ ^1 

(k+1) 

(k) 

Jk) 

®m+l 

°‘q-k-3^m 

^m+1 * * ’ ■ ’ ^ 

(k+1) 



*k+2 

" ^q-k-3 ^k+1 

(q-2) 

where k=0, 1 , . , 

. q-3. The desired coefficients are then given by ^q, a^ , 


m=l, 2, ..., q-1. The operations in (5) require only a single vector of 

(k) 

storage since and . k=0, 1, ... can all occupy the same location. 

Denoting this location by we nave 

'^q-1 " “q-2'^q-l 


'^‘q+m-k-2 '*q-k-3'^q+m-k-2 "*^'^'q+m-k-l , m=0,l, ..., k 

’^q-l 'q-k-3'^q-l 



1 . 


. q-3 
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One gets a similar recurrence from eq. (5.11) of [3] by replacing the itera- 
tion indices k=l, 2, .... i; i=q-l, .... 1 with i=q-l, k; k=l. 2, .... q-I, 
then replacing i with q+k-j-1 and then interchanging the order of iteration 
indices once again. 

For the case of converting to the Chebyshev basis, we have A(m) = 0, 

B(0) = -1, B(m) = -2, m > 0, C(m) = 1. Substituting as before, and then setting 

.(k) 


<4 /Vk-2> = 4^’- ■ aI'*. .. > 0. 


m '“q-k-2' ■ ’n 


one obtains 


(0) 


II 

o 

<ro 

1 -^ 

^.2^.1 

(k+1) 

(k) (k) 

"o 

“ “q-k-3 [2*q.k-2 ^4 * h ' '^q-3*q. 

(k+1) 

Jk) Jk) (k) 


° (’■'^>“q-k-3 l-Vl * Vl^ * 2 

(k+1) 

Jk) (k) 

®k 

- <’/2)“q-k-3 Vi * 4 

Jk+1) 

(k) 

®k+l 

= <’/2)“q-k-3 4 

where k=0, 1 , . , 

{' 

. q-3, and the desired coefficients are a 

m 

(q-1) 

(q-2) 

o 

«o 

CO 

= ♦o + (1/2) a, 

(q-1) 

Jq-2) Jq-2) 

^m 

' (1/2) ], m=l, 2 q-3 

^(q-1) 

Jq-2) 

^-2 

= (1/2) 3^.3 

(q-1 ) 

(q-2) 

4-1 

= (1/2) 4-2 


^( 0 ) 

‘o 


(q-I) 


given by 
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As before, the cdiculations in (7) can be done in a single vector 
of storage as follows (initially 

(5) = 2aq_2ii^q.-, 

^q-2 " ^°‘q-3'^q-2 '^q-1 

•^q-1 = 

’^q-k-2 "" “‘q-k-a ^^'^q-k-2 '^q-k^ '^q-k-1 

’^q-k-2+m ' ^^'*'^^^q-k-3^'^q-k-2+m ^ '^q-k+m^ '^q-k-l+m 

m=l , 2, .... k-1 

^ k-1, 2, .... q-3 

♦q-2 ' <’''^>“q-k-3^-2 * *q-l ! 

♦q-1 ° (’/^'“q-k-a^q-l 

The final Chebyshev coefficients are then obtained, using (8) as follows 
(10) 'i'o ^ (V2)<l^-j 

tq.2 = 0/2)4^q,2 
'^q-1 = 

If one is using the method outlined in [3] or the code in [4], the 
a's are available at the time of integration. Clearly some time can be saved 
in the above calculations if t is redefined to be (t-t^)/2h, since this has 
the effect of setting = -2ot^ and as a result the multiplications by 1/2 are 
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eliminated. {For this case, simply replace a. by in (9), and replace all 
1/2's in (10) with 1. For t = ^ everyplace we just indicated 

a multiplication by 2.) 

If one wants to approximate the solution of a differential equation 
by a Chebyshev polynomial, and thus wants to get the coefficients for the in- 
tegrated polynomial, one can save a 1 ittle computation by computing the integra- 
ted coefficients after the calculations indicated in (7) or (9). Thus for 
3 « n ^ q-4 one has for 3^, the n-th coefficient o^ the integrated polynomial 


01 ) 


3. 


(q-n 

(q-1) 


(l/2n)(a„., 

■ ^n+1 

) after (8) 

(q-2) 

(q-2) 


n/4n)(a„.2 

■ ^n+2 

) after (7) 



Computing Coefficients for the Integrated or Differentiated Polynomial 


When integrating stiff equations, the best integration method is 
different depending on how far the most active component is from its local 
equilibrium point. If a procedure such as that in [3] is used, changes in 
the method involve either getting coefficients for the integrated polynomial 
when leaving a transient region, or for the differentiated polynomial when a 
transient is encountered. We have also encountered a desire for such a 
transformation in another context. 


The best general algorithm we have been able to find for the case 
of variable stepsize involves running through the recurrences in (6), repla- 
cing with k=q, ..., 1 for integration (or with (k+l , k-0, 1, 

.... q-2 for differentiation). For integration, is set to the constant of 
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integration. Then simply run th3 recurrence (6) in reverse order 

( 12 ) 

’^q'+m-k-2 ^ ^'^q'+m-k-l " 'l'q'+m-k-2^^“q'-k-3» ° j k=q'-2. 

^•-1 = ^'-l/“q'.k-3 

where q' = q+1 , for integration, and q'=q-l for differentiation. In the case 
of a constant stepsize, one can accomplish this goal more efficiently by re- 
peated differencing of the corrector equation 


q-1 ^ 

(13) Vy = h E YfV y' 
k=0 


Thus 


q“>i k+j-i 

(14) V'^y^ = h YjV y'^ , j = l, 2, .... q. 

k=n 


and the reverse operation simply involves treating eq. (14) as a triangular 

k 

system to be solved for the V y'^. 
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